home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Pascal Super Library
/
Pascal Super Library (CW International)(1997).bin
/
MATH
/
MATH2
/
BETAINV.PAS
< prev
next >
Wrap
Pascal/Delphi Source File
|
1985-05-17
|
9KB
|
261 lines
(*--------------------------------------------------------------------------*)
(* BetaInv -- Inverse Beta Distribution *)
(*--------------------------------------------------------------------------*)
FUNCTION BetaInv( P, Alpha, Beta : REAL;
MaxIter: INTEGER;
Dprec: INTEGER;
VAR Iter: INTEGER;
VAR Cprec: REAL;
VAR Ierr: INTEGER ) : REAL;
(*--------------------------------------------------------------------------*)
(* *)
(* Function: BetaInv *)
(* *)
(* Purpose: Calculates inverse Beta distribution *)
(* *)
(* Calling Sequence: *)
(* *)
(* B := BetaInv( P, Alpha, Beta : REAL *)
(* MaxIter : INTEGER; *)
(* Dprec : INTEGER; *)
(* VAR Iter : INTEGER; *)
(* VAR Cprec : REAL; *)
(* VAR Ierr : INTEGER ) : REAL; *)
(* *)
(* P --- Cumulative probability for which percentage *)
(* point is to be found. *)
(* Alpha --- First parameter of Beta distribution *)
(* Beta --- Second parameter of Beta distribution *)
(* MaxIter --- Maximum iterations allowed *)
(* Dprec --- no. of digits of precision desired *)
(* Iter --- no. of iterations actually performed *)
(* Cprec --- no. of digits of precision actually obtained *)
(* Ierr --- error indicator *)
(* = 0: OK *)
(* = 1: argument error *)
(* = 2: evaluation error during iteration *)
(* *)
(* B --- Resulting percentage point of Beta dist. *)
(* *)
(* Calls: CDBeta (cumulative Beta distribution) *)
(* *)
(* Remarks: *)
(* *)
(* Because of the relationship between the Beta *)
(* distribution and the F distribution, this *)
(* routine can be used to find the inverse of *)
(* the F distribution. *)
(* *)
(* The percentage point is returned as -1.0 *)
(* if any error occurs. *)
(* *)
(* Newtons' method is used to search for the percentage point. *)
(* *)
(* The algorithm is based upon AS110 from Applied Statistics. *)
(* *)
(*--------------------------------------------------------------------------*)
VAR
Eps : REAL;
Xim1 : REAL;
Xi : REAL;
Xip1 : REAL;
Fim1 : REAL;
Fi : REAL;
W : REAL;
Cmplbt : REAL;
Adj : REAL;
Sq : REAL;
R : REAL;
S : REAL;
T : REAL;
G : REAL;
A : REAL;
B : REAL;
PP : REAL;
H : REAL;
A1 : REAL;
B1 : REAL;
Eprec : REAL;
Done : BOOLEAN;
Jter : INTEGER;
LABEL 10, 30, 9000;
BEGIN (* BetaInv *)
Ierr := 1;
BetaInv := P;
(* Check validity of arguments *)
IF( ( Alpha <= 0.0 ) OR ( Beta <= 0.0 ) ) THEN GOTO 9000;
IF( ( P > 1.0 ) OR ( P < 0.0 ) ) THEN GOTO 9000;
(* Check for P = 0 or 1 *)
IF( ( P = 0.0 ) OR ( P = 1.0 ) ) THEN
BEGIN
Iter := 0;
Cprec := MaxPrec;
GOTO 9000;
END;
(* Set precision *)
IF Dprec > MaxPrec THEN
Dprec := MaxPrec
ELSE IF Dprec <= 0 THEN
Dprec := 1;
Cprec := Dprec;
Eps := PowTen( -2 * Dprec );
(* Flip params if needed so that *)
(* P for evaluation is <= .5 *)
IF( P > 0.5 ) THEN
BEGIN
A := Beta;
B := Alpha;
PP := 1.0 - P;
END
ELSE
BEGIN
A := Alpha;
B := Beta;
PP := P;
END;
(* Generate initial approximation. *)
(* Several different ones used, *)
(* depending upon parameter values. *)
Ierr := 0;
Cmplbt := ALGama( A ) + ALGama( B ) - ALGama( A + B );
Fi := Ninv( 1.0 - PP );
IF( ( A > 1.0 ) AND ( B > 1.0 ) ) THEN
BEGIN
R := ( Fi * Fi - 3.0 ) / 6.0;
S := 1.0 / ( A + A - 1.0 );
T := 1.0 / ( B + B - 1.0 );
H := 2.0 / ( S + T );
W := Fi * SQRT( H + R ) / H - ( T - S ) *
( R + 5.0 / 6.0 - 2.0 / ( 3.0 * H ) );
Xi := A / ( A + B * EXP( W + W ) );
END
ELSE
BEGIN
R := B + B;
T := 1.0 / ( 9.0 * B );
T := R * PowerI( ( 1.0 - T + Fi * SQRT( T ) ) , 3 );
IF( T <= 0.0 ) THEN
Xi := 1.0 - EXP( ( LN( ( 1.0 - PP ) * B ) + Cmplbt ) / B )
ELSE
BEGIN
T := ( 4.0 * A + R - 2.0 ) / T;
IF( T <= 1.0 ) THEN
Xi := EXP( (LN( PP * A ) + Cmplbt) / PP )
ELSE
Xi := 1.0 - 2.0 / ( T + 1.0 );
END;
END;
(* Force initial estimate to *)
(* reasonable range. *)
IF ( Xi < 0.0001 ) THEN Xi := 0.0001;
IF ( Xi > 0.9999 ) THEN Xi := 0.9999;
(* Set up Newton-Raphson loop *)
A1 := 1.0 - A;
B1 := 1.0 - B;
Fim1 := 0.0;
Sq := 1.0;
Xim1 := 1.0;
Iter := 0;
Done := FALSE;
(* Begin Newton-Raphson loop *)
REPEAT
Iter := Iter + 1;
Done := Done OR ( Iter > MaxIter );
Fi := CDBeta( Xi, A, B, Dprec+1, MaxIter, Eprec, Jter, Ierr );
IF( Ierr <> 0 ) THEN
BEGIN
Ierr := 2;
Done := TRUE;
END
ELSE
BEGIN
Fi := ( Fi - PP ) * EXP( Cmplbt + A1 * LN( Xi ) +
B1 * LN( 1.0 - Xi ) );
IF( ( Fi * Fim1 ) <= 0.0 ) THEN Xim1 := Sq;
G := 1.0;
10: REPEAT
Adj := G * Fi;
Sq := Adj * Adj;
IF( Sq >= Xim1 ) THEN G := G / 3.0;
UNTIL( Sq < Xim1 );
Xip1 := Xi - Adj;
IF( ( Xip1 < 0.0 ) OR ( Xip1 > 1.0 ) ) THEN
BEGIN
G := G / 3.0;
GOTO 10;
END;
IF( Xim1 <= Eps ) THEN GOTO 30;
IF( Fi * Fi <= Eps ) THEN GOTO 30;
IF( ( Xip1 = 0.0 ) OR ( Xip1 = 1.0 ) ) THEN
BEGIN
G := G / 3.0;
GOTO 10;
END;
IF( Xip1 <> Xi ) THEN
BEGIN
Xi := Xip1;
Fim1 := Fi;
END
ELSE
Done := TRUE;
End;
UNTIL( Done );
30:
BetaInv := Xi;
IF( P > 0.5 ) THEN BetaInv := 1.0 - Xi;
IF ABS( Xi - Xim1 ) <> 0.0 THEN
Cprec := -LogTen( ABS( Xi - Xim1 ) )
ELSE
Cprec := MaxPrec;
9000:
IF Ierr <> 0 THEN BetaInv := -1.0;
END (* BetaInv *);